Mouse small intestine analysis, 20 classes
Jason T. Serviss
2019-08-22
- add table of connections found per number of cells in multiplet
- Fig 1: Classes
- Fig 2: Cell type gene expression
- Fig 3: Cell cycle
- Fig 4: Connections per multiplet
- Fig 5: Fraction histogram
- Fig 6: Detected cell types vs. cost
- Fig 7: Estimated cell numbers vs. cost
- Fig 8: Estimated cell number vs. Detected cell number
- Fig 9: Detected cell number vs. Total counts
- Fig 10: Detected cell number vs. Total ERCC counts
- Fig 11: Relative frequency of cell types in singlets vs. deconvoluted multiplets
- Fig 12: All connections
- Fig 13: Filtered
add table of connections found per number of cells in multiplet
packages <- c("CIMseq", "CIMseq.testing", "tidyverse", "circlize", "printr")
purrr::walk(packages, library, character.only = TRUE)
rm(packages)
##DATA
load('../data/CIMseqData.rda')
load('../data/sObj.rda')
if(!dir.exists('../figures')) dir.create('../figures')
#there are 2 cells that were classified as colon but sorted as SI. These have to
#be removed manually
c <- getData(cObjSng, "classification")
s <- names(c[c %in% c("8", "13")])
i <- which(colnames(getData(cObjSng, "counts")) %in% s)
cObjSng <- CIMseqSinglets(
getData(cObjSng, "counts")[, -i],
getData(cObjSng, "counts.ercc")[, -i],
getData(cObjSng, "dim.red")[-i, ],
getData(cObjSng, "classification")[-i]
)Fig 1: Classes
p <- plotUnsupervisedClass(cObjSng, cObjMul, palette('si'))
p
ggsave(
plot = p,
filename = '../figures/MGA.enge20SI.classes.pdf',
device = cairo_pdf,
height = 240,
width = 240,
units = "mm"
)Fig 2: Cell type gene expression
p <- plotUnsupervisedMarkers(
cObjSng, cObjMul,
c("Lgr5", "Muc2", "Ptprc", "Chga", "Alpi", "Lyz1", "Dclk1"),
pal = RColorBrewer::brewer.pal(8, "Set1")
)
p
ggsave(
plot = p,
filename = '../figures/MGA.enge20SI.markers.pdf',
device = cairo_pdf,
height = 240,
width = 240,
units = "mm"
)Fig 3: Cell cycle
p <- plotUnsupervisedMarkers(
cObjSng, cObjMul, c("Mki67"),
pal = RColorBrewer::brewer.pal(8, "Set1")
)
p
ggsave(
plot = p,
filename = '../figures/MGA.enge20SI.cellcycle.pdf',
device = cairo_pdf,
height = 240,
width = 240,
units = "mm"
)Fig 4: Connections per multiplet
adj <- adjustFractions(cObjSng, cObjMul, sObj)
table(apply(adj, 1, sum))| 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| 19 | 80 | 217 | 86 | 31 | 2 |
Fig 5: Fraction histogram
tibble(fractions = c(fractions)) %>%
ggplot() +
geom_histogram(aes(fractions), binwidth = 0.01) +
theme_bw()
Range of fractions picked after adjustment.
range(fractions[adj == 1])## [1] 0.02485269 0.99998969Fig 6: Detected cell types vs. cost
tibble(
nCellTypes = apply(adj, 1, sum),
cost = getData(sObj, "costs")
) %>%
ggplot() +
geom_boxplot(aes(nCellTypes, cost, group = nCellTypes)) +
scale_x_continuous(name = "Detected cell types", breaks = 0:max(apply(adj, 1, sum))) +
theme_bw()
Fig 7: Estimated cell numbers vs. cost
tibble(
sample = rownames(getData(sObj, "fractions")),
cost = unname(getData(sObj, "costs"))
) %>%
inner_join(
select(estimateCells(cObjSng, cObjMul), sample, estimatedCellNumber),
by = "sample"
) %>%
mutate(estimatedCellNumber = round(estimatedCellNumber)) %>%
ggplot() +
geom_boxplot(aes(estimatedCellNumber, cost, group = estimatedCellNumber)) +
scale_x_continuous(
name = "ERCC estimated cell number",
breaks = 0:max(round(pull(estimateCells(cObjSng, cObjMul), estimatedCellNumber)))
) +
theme_bw()
Fig 8: Estimated cell number vs. Detected cell number
ercc <- filter(estimateCells(cObjSng, cObjMul), sampleType == "Multiplet")
nConnections <- apply(adj, 1, sum)
nConnections <- nConnections[match(ercc$sample, names(nConnections))]
tibble(
detectedConnections = round(nConnections),
estimatedCellNumber = round(ercc$estimatedCellNumber)
) %>%
ggplot() +
geom_boxplot(aes(estimatedCellNumber, detectedConnections, group = estimatedCellNumber)) +
scale_x_continuous(
name = "ERCC estimated cell number",
breaks = 0:max(round(ercc$estimatedCellNumber))
) +
scale_y_continuous(
name = "Detected cell number",
breaks = 0:max(round(nConnections))
) +
theme_bw()
Fig 9: Detected cell number vs. Total counts
tibble(
sample = names(nConnections),
detectedConnections = nConnections
) %>%
inner_join(tibble(
sample = colnames(getData(cObjMul, "counts")),
total.counts = colSums(getData(cObjMul, "counts"))
), by = "sample") %>%
ggplot() +
geom_boxplot(aes(detectedConnections, total.counts, group = detectedConnections)) +
scale_x_continuous(
name = "Detected cell number",
breaks = 0:max(nConnections)
) +
scale_y_continuous(name = "Total counts") +
theme_bw()
Fig 10: Detected cell number vs. Total ERCC counts
tibble(
sample = names(nConnections),
detectedConnections = nConnections
) %>%
inner_join(tibble(
sample = colnames(getData(cObjMul, "counts")),
total.ercc = colSums(getData(cObjMul, "counts.ercc"))
), by = "sample") %>%
ggplot() +
geom_boxplot(aes(detectedConnections, total.ercc, group = detectedConnections)) +
scale_x_continuous(
name = "Detected cell number",
breaks = 0:max(nConnections)
) +
scale_y_continuous(name = "Total ERCC counts") +
theme_bw()
Fig 11: Relative frequency of cell types in singlets vs. deconvoluted multiplets
singlets <- c(table(getData(cObjSng, "classification")))
singlets <- singlets / sum(singlets)
deconv <- colSums(adjustFractions(cObjSng, cObjMul, sObj))
deconv <- deconv[match(names(singlets), names(deconv))]
deconv <- deconv / sum(deconv)
if(!identical(names(singlets), names(deconv))) stop("name mismatch")
p <- tibble(
class = names(singlets),
singlet.freq = singlets,
multiplet.freq = deconv
) %>%
ggplot() +
geom_point(aes(singlet.freq, multiplet.freq, colour = class), size = 3) +
scale_colour_manual(values = palette('si')[order(names(palette('si')))]) +
xlim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
ylim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
geom_abline(slope = 1, intercept = 0, lty = 3, colour = "grey") +
labs(x = "Singlet relative frequency", y = "Multiplet relative frequency") +
guides(colour = guide_legend(title = "Cell Type")) +
theme_bw()
p
ggsave(
plot = p,
filename = '../figures/MGA.enge20SI.sngMulRelFreq.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)Fig 12: All connections
plotSwarmCircos(
sObj, cObjSng, cObjMul, classOrder = classOrder.MGA('si'), classColour = palette('si')[classOrder.MGA('si')],
h.ratio = 0.85
)## Joining, by = "class"
Fig 13: Filtered
Only detected duplicates, triplicates, and quadruplicates.
ERCC estimated cell number set to max 4.
Weight cutoff = 10.
# adj <- adjustFractions(cObjSng, cObjMul, sObj, binary = TRUE)
# samples <- rownames(adj)
# rs <- rowSums(adj)
# keep <- rs == 2 | rs == 3 | rs == 4
plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 10,
classOrder = classOrder.MGA('si'), theoretical.max = 4, classColour = palette('si')[classOrder.MGA('si')],
h.ratio = 0.85, alpha = 1e-3
)## Joining, by = "class"
pdf('../figures/MGA.enge20SI.circos.pdf', width = 9.5, height = 9.5, onefile = FALSE)
plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 10,
classOrder = classOrder.MGA('si'), theoretical.max = 4, classColour = palette('si')[classOrder.MGA('si')],
h.ratio = 0.85, alpha = 1e-3
)## Joining, by = "class"dev.off()## quartz_off_screen
## 2Calculate probablity of paneth - other cell type interaction as the fraction of other cell types observed in multiplets reported to contain a paneth cell.
pdata <- adjustFractions(cObjSng, cObjMul, sObj, theoretical.max = 4) %>%
matrix_to_tibble("sample") %>%
filter(Paneth == 1) %>%
select(-Paneth) %>%
gather(class, binary, -sample) %>%
group_by(sample) %>%
summarize(others = paste(class[binary == 1], collapse = ", ")) %>%
mutate(others = map(others, ~str_split(.x, ", ")[[1]])) %>%
unnest() %>%
filter(others != "") %>%
group_by(others) %>%
summarize(prob = n() / nrow(.)) %>%
rename(class = others) %>%
full_join(tibble(class = unique(getData(cObjSng, "classification")))) %>%
filter(class != "Paneth") %>%
replace_na(list(prob = 0))## Joining, by = "class"p <- pdata %>%
ggplot() +
geom_bar(aes(class, prob), stat = "identity", position = position_dodge(width = 1)) +
geom_text(aes(class, prob + 0.01, label = round(prob, digits = 3))) +
theme_bw() +
labs(y = "Probability") +
theme(axis.title.x = element_blank())
p
ggsave(
plot = p,
filename = '../figures/MGA.enge20.PanethIntProb.pdf',
device = cairo_pdf,
height = 240,
width = 240,
units = "mm"
)Calculate the probability of observing Lgr5 expression in multiplets that express Lyz1.
#calculate cutoff for Lyz1 based on singlets
cut <- getData(cObjSng, "counts.cpm") %>%
.['Lyz1', ] %>%
tibble(sample = names(.), Lyz1 = .) %>%
filter(getData(cObjSng, "classification") == "Paneth") %>%
pull(Lyz1) %>%
min()
# p <- getData(cObjMul, "counts.cpm") %>%
# .[c("Lyz1", "Lgr5"), ] %>%
# t() %>%
# matrix_to_tibble("sample") %>%
# #filter(Lyz1 > cut) %>% #include only Lyz1 positive
# mutate(
# express.lgr5 = if_else(Lgr5 > 0, TRUE, FALSE),
# express.lyz1 = if_else(Lyz1 > cut, TRUE, FALSE)
# ) %>%
# count(express.lgr5, express.lyz1) %>%
# group_by(express.lyz1) %>%
# mutate(total = sum(n)) %>%
# mutate(lgr5.prob = n / total) %>%
# ungroup() %>%
# filter(express.lgr5) %>%
# ggplot() +
# geom_bar(aes(express.lyz1, lgr5.prob), stat = "identity", position = position_dodge(width = 1)) +
# labs(x = "Lyz1 expressed", y = "Lgr5 probability") +
# ggthemes::theme_few()
getData(cObjMul, "counts.cpm") %>%
.[c("Lyz1", "Lgr5"), ] %>%
t() %>%
matrix_to_tibble("sample") %>%
filter(Lyz1 > cut) %>%
mutate(express.lgr5 = if_else(Lgr5 > 0, TRUE, FALSE)) %>%
count(express.lgr5) %>%
mutate(total = sum(n)) %>%
filter(express.lgr5) %>%
mutate(prob = n / total) %>%
pull(prob)## [1] 0.9sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] printr_0.1 circlize_0.4.6 forcats_0.4.0
## [4] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2
## [7] readr_1.3.1 tidyr_0.8.3 tibble_2.1.3
## [10] ggplot2_3.2.1 tidyverse_1.2.1 CIMseq.testing_0.0.2
## [13] CIMseq_0.3.0.2
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-141 matrixStats_0.54.0 lubridate_1.7.4
## [4] RColorBrewer_1.1-2 gmodels_2.18.1 httr_1.4.1
## [7] tools_3.6.1 backports_1.1.4 R6_2.4.0
## [10] lazyeval_0.2.2 BiocGenerics_0.30.0 colorspace_1.4-1
## [13] withr_2.1.2 tidyselect_0.2.5 gridExtra_2.3
## [16] compiler_3.6.1 cli_1.1.0 rvest_0.3.4
## [19] xml2_1.2.1 labeling_0.3 scales_1.0.0
## [22] digest_0.6.20 rmarkdown_1.14 pkgconfig_2.0.2
## [25] htmltools_0.3.6 highr_0.8 rlang_0.4.0
## [28] GlobalOptions_0.1.0 ggthemes_4.2.0 readxl_1.3.1
## [31] rstudioapi_0.10 shape_1.4.4 farver_1.1.0
## [34] generics_0.0.2 jsonlite_1.6 gtools_3.8.1
## [37] magrittr_1.5 Rcpp_1.0.2 munsell_0.5.0
## [40] S4Vectors_0.22.0 viridis_0.5.1 stringi_1.4.3
## [43] yaml_2.2.0 ggraph_1.0.2 MASS_7.3-51.4
## [46] Rtsne_0.15 plyr_1.8.4 grid_3.6.1
## [49] parallel_3.6.1 gdata_2.18.0 listenv_0.7.0
## [52] ggrepel_0.8.1 crayon_1.3.4 lattice_0.20-38
## [55] haven_2.1.1 hms_0.5.0 zeallot_0.1.0
## [58] knitr_1.23 pillar_1.4.2 igraph_1.2.4.1
## [61] pso_1.0.3 future.apply_1.3.0 codetools_0.2-16
## [64] stats4_3.6.1 glue_1.3.1 evaluate_0.14
## [67] modelr_0.1.4 vctrs_0.2.0 tweenr_1.0.1
## [70] cellranger_1.1.0 gtable_0.3.0 RANN_2.6.1
## [73] polyclip_1.10-0 future_1.14.0 assertthat_0.2.1
## [76] xfun_0.8 gridBase_0.4-7 ggforce_0.2.2
## [79] broom_0.5.2 tidygraph_1.1.2 viridisLite_0.3.0
## [82] globals_0.12.4